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ABSTRACT 



o 

in 
c3 



We study the nonlinear evolution of the Rossby wave instability in thin disks 



o 
o, 

using global 2D hydrodynamic simulations. The detailed linear theory of this 
nonaxisymmetric instability was developed earlier by Lovelace et al. and Li et al., 
' who found that the instability can be excited when there is an extremum in the radial 

Q ■ profile of an entropy-modified version of potential vorticity. The key questions we 

£Sj ■ are addressing in this paper are: (1) What happens when the instability becomes 

nonlinear? Specifically, does it lead to vortex formation? (2) What is the detailed 
behavior of a vortex? (3) Can the instability sustain itself and can the vortex last 
a long time? Among various initial equilibria that we have examined, we generally 
\ find that there are three stages of the disk evolution: (1) The exponential growth of 

| the initial small amplitude perturbations. This is in excellent agreement with the 

linear theory; (2) The production of large scale vortices and their interactions with 
the background flow, including shocks. Significant accretion is observed due to these 
■ vortices. (3) The coupling of Rossby waves /vortices with global spiral waves, which 

facilitates further accretion throughout the whole disk. Even after more than 20 
revolutions at the radius of vortices, we find that the disk maintains a state that is 
populated with vortices, shocks, spiral waves/shocks, all of which transport angular 
momentum outwards. We elucidate the physics at each stage and show that there 
is an efficient outward angular momentum transport in stages (2) and (3) over most 
parts of the disk, with an equivalent Shakura-Sunyaev angular momentum transport 
parameter a in the range from 10 -4 to 10~ 2 . By carefully analyzing the flow structure 
around a vortex, we show why such vortices prove to be almost ideal "units" in 
transporting angular momentum outwards, namely by positively correlating the radial 
and azimuthal velocity components. In converting the gravitational energy to the 
internal energy, we find some special cases in which entropy can remain the same while 
angular momentum is transported. This is different from the classical a disk model 
which results in the maximum dissipation (or entropy production). The dependence 
of the transport efficiency on various physical parameters are examined and effects of 
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radiative cooling are briefly discussed as well. We conclude that Rossby wave/vortex 
instability is an efficient, purely hydrodynamic mechanism for angular momentum 
transport in thin disks, and may find important applications in many astrophysical 
systems. 



Subject headings: Accretion Disks — Hydrodynamics — Instabilities — Waves 



1. Introduction 

Understanding the physics of accretion disks has remained a great challenge in astrophysics 
for decades. Matter has to lose angular momentum in order to fall deeper into a gravitational 
potential well. The release of the gravitational binding energy then becomes one of the most 
powerful energy sources in the universe. Various models for angular momentum transport have 
been proposed, including those having a purely radial transport (i.e., within the disk) and those 
using outflows (e.g., MHD jets). One promising mechanism for removing angular momentum 
locally within the disk is via MHD turbulence in disks ( [Balbus fe Hawley 199S| ). Disks must be 
made of relatively hot and sufficiently ionized plasma for this mechanism to operate, however, 
because a strong coupling between magnetic fields and plasma is required. On the other hand, 
there are several types of astrophysical disks where such conditions are not fully satisfied. Thus, a 
purely hydrodynamic means of angular momentum transport is still needed (see papaloizou fc Lin 



1995 for a recent review) 



In two previous papers, Lovelace et al. (2000, hereafter Paper I) and Li et al. (2000, hereafter 
Paper II) have presented a detailed analysis of the linear theory of a global, nonaxisymmetric 
hydrodynamic instability in thin (2D) disks. The disk becomes unstable when the conditions of 
Rayleigh's inflection point theorem are violated, which is indicated by the radial profile of a key 
function C(r) = (SO/ k?)S 2 ^ . This function is an entropy-modified version of potential vorticity. 
Here, S(r) is the surface mass density of the disk, O(r) the angular rotation rate, S(r) the specific 
entropy, T the adiabatic index, and K,(r) the radial epicyclic frequency. It has been shown that 
a sufficient variation of pressure over a length scale that is a few times the disk thickness can 
cause the disk to be unstable to nonaxisymmetric perturbations even though the disk is still 
locally stable to axisymmetric perturbations according to the so-called Rayleigh determinant 
or the Solberg-Hoiland criterion when the pressure effects are included (Paper II). The linear 
theory shows that the unstable modes have a dispersion relation similar to that of Rossby waves 
in atmospheric studies (Paper I). The term "Rossby wave instability" (RWI) was introduced in 
that paper. The dependence of RWI on various physical parameters has been examined in Paper 
II and its relations with other hydrodynamic instabilities, especially the Papaloizou & Pringle 



instability ( Papaloizou fc Pringle 1985 ), have been discussed in Paper II as well. More generally, 
since the pioneering work by Papaloizou & Pringle (1985), nonaxisymmetric instabilities in disks 
have received an enormous amount of attention with various degrees of success. In particular, the 
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important role of "vortensity" (vorticity divided by surface density) in determining the stability 
of the disk was first discussed in Lovelace & Holfeld (1978) and was studied in greater detail 
in Papaloizou & Lin (1989). Nonaxisymmetric convective (in the vertical direction) instability 
is also explored in Lin, Papaloizou, & Kley (1993). Sellwood & Kahn (1991) demonstrated that 
a "narrow groove" in the angular momentum density profile can be very unstable. Interestingly, 
Toomre (1981) had already noticed certain unstable modes associated with disk edges (termed 
"edge modes"). Even though these earlier studies have mostly used a homentropic equation of 
state (i.e., effects due to an entropy gradient were usually not present), they are, in general, 
consistent with the above mentioned criterion by having an inflection point in the radial profile of 
the key function C{r). 

There have been several other recent studies of the role of Rossby waves in accretion disks. 
Sheehan et al. (1999) performed a linear analysis of the generation of Rossby waves and their 
propagation in protoplanetary disks. They furthermore speculated on the possibility of forming 
vortices and zonal jets (similar to planetary atmosphere dynamics). Using extensive nonlinear 
disk simulations, Godon & Livio (2000) have investigated the formation of vortices in a viscous, 
compressible flow and Klahr & Bodenheimer (2000) have shown that vortices can be produced by a 
global radial entropy gradient, possibly via a baroclinic instability from which angular momentum 
is transported outwards with an efficiency at a ~ 10~ 3 level. 

Besides the possible role of Rossby waves/ vortices, another important purely hydrodynamic 
angular momentum transport mechanism is through spiral waves/shocks. This has been proposed 



for systems such as cataclysmic variables (CVs; see Spruit 199l| for a review) and accreting neutron 
stars ( [Michel 1984| ), where the nonaxisymmetry in the disk flow is caused by external torques 
acting on the disk either from a close companion or an asymmetric central rotating body (such as 
the magnetosphere of a neutron star). In the linear theory analysis of Tagger & Pellat (1999), a 
possible mechanism of coupling Rossby waves with spiral (density) waves through the corotation 
resonance was briefly discussed. 

In this paper, we use extensive 2D hydrodynamic simulations to study the nonlinear evolution 
of RWI, based on the knowledge we have gained from the linear theory analysis. After a brief 
description of our 2D hydro code, we show how the initial states of disk simulations are set up 
in §^, We then present the simulation results in §|| and §|j. In §|5], we discuss several important 
physical issues that are associated with this work. Conclusions are given in §||. 



2. The 2D Hydrodynamic Model 

Several simplifying assumptions are employed in this study. The disk is assumed to be 
geometrically thin so that the hydro equations can be reduced to 2D with vertically integrated 
quantities. The effects due to magnetic field and self-gravity are omitted and the Newtonian 
potential is used throughout our simulations. The disk is treated as an isolated system so that 
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there is no externally supplied mass inflow, but mass outflow through the disk radial boundaries 
(both sides) is allowed. 



2.1. The differential equations 



The Euler equations are the governing equations for our 2D (in r, (f> plane) , inviscid disk flow 
with a central gravitating object. The usual variables in Euler equations are u = {£, v r , v^, E}, 
where v r ,v,p are the radial and azimuthal velocities, respectively, E is the total energy 
E = P/(r — 1) + 0.5 x £(f 2 + vfy, T is the ideal gas adiabatic index, and P is the vertically 
integrated pressure. Here, we choose to use a new set of variables u = {rS, rT,v r , r 2 Su^, rE}. 
One advantage of this choice is to eliminate the nonzero source term in the angular momentum 
equation. Consequently, the Euler equations in cylindrical coordinates become 



where 
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The S/r term in S is the normalized central gravitational potential. The zero components of S 
express the conservation of mass and angular momentum. 
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2.2. The Initial Conditions 



We study the evolution of RWI by first setting up the initial equilibria which are steady, 
axisymmetric (d/dcj) = 0) and with zero radial velocity (v r = 0). Following the analysis in Paper 
II, we first specify the surface density S(r) and temperature T(r) if the disk does not have a 
constant entropy initially As the main focus of this paper is to study the nonlinear evolution 
of the linear instability found in Paper I and II, we follow the initial setups in Paper II and 
concentrate on two types of initial equilibrium. One has a Gaussian-shaped "bump" in pressure 
and the other has a step-jump in pressure. For the sake of completeness, we recap the functions 
used to describe the bump/jump (Paper II), 



Bg(r) = 1 + (A- l)exp 



for the Gaussian bump and 



Bj{r) 



1 + 4 

2 



tanh 



r — ro 
Ar 



Ar 



+ 1 



(6) 



(7) 



for the step jump. Quantities A and Ar measure the amplitude and width of the bump /jump 
respectively, and ro is radial location of the bump/jump. Specifically, we have considered 4 basic 
types, 




= Po Bg 

= P [p(r)/po] T ' 

Po (r/r )~ 3/4 
T (r/r )- 3 / 4 Bg 
PoTo (r/r )- 3 / 2 Bg 



(8) 



(9) 



HSJ : 



p(r) 



Po Bj 



P(r) = P [p(r)/po} 1 



(10) 



p( r ) = Po {r/r ) 3/4 Bj 
NS.T : I T(r) = T (r/r )- 3 / 4 Bj , (11) 

P(r) = p T (r/r )- 3 / 2 Bj 2 

which are named the homentropic Gaussian bump (HGB), nonhomentropic Gaussian bump 
(NGB), homentropic step jump (HSJ), and nonhomentropic step jump (NSJ) cases, respectively. 
These equations have an overall normalization so that c 2 = TTo/So = TTo = 0.01v 2 ifc (ro). Note 
that even though we have used either P oc ST or P oc S r in obtaining the initial pressure 
distributions, we only use P oc ST as the equation of state during the subsequent disk evolution. 
For a given S(r) and P(r), we use the radial force balance to calculate the azimuthal velocity v^, 
which is very close to Keplerian velocities except the slight modification by the pressure gradient 
(Paper II). 
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To be consistent with the 2D approximation, we require that the length scale of the pressure 
variation L p = TP /\dP / dr\ (in units of tq) is larger than the disk scale height ~ Cs/v^ at r$. Note 
that L p is always larger than Ar. For example, using equation ([?]), if Ar/r = 0.05 and A = 0.65, 
the minimum of L p /tq is ~ 0.2, which is twice the thickness of disk when c s jv^ at tq is 0.1. In the 
following runs, we have used both Ar/ro = 0.05 and 0.1. The amplitude A is the main parameter 
to be varied. As discussed in Paper II, the RWI growth rate 7 of RWI is a strong function of c s 
and L p , and the growth rate 7 becomes a fraction (~ 0.2) of the Keplerian rotation O(ro) when 
L p is ~ 2x the disk thickness. This is because the pressure gradient is the only force available to 
perturb the Keplerian flow in an inviscid disk. 

There are several considerations that go into constructing these different types. The disks of 
HGB and HSJ types start with a single entropy for the whole disk. They are relatively "clean" 
systems, and so help us to determine whether an entropy gradient is needed in the development 
of RWI. The linear theory in Paper II says dS/dr can be zero for RWI to grow. We can test 
it and compare them with the cases of NGB and NSJ where a background radial gradient in 
entropy is present. The power law dependence of S(r) and T(r) in NGB and NSJ mimics the 
distribution from a steady state a— disk model, although the code can handle an arbitrary slope. 
Furthermore, the bump in NGB is in temperature only whereas NSJ has a jump in both density 
and temperature. 

2.3. Description of the numerical method 

The system (|]) is integrated using a dimensional-splitting method. Furthermore, a local 
co- moving frame for the <f> sweeps is employed. It was observed by Masset (2000) that this reduces 
the computing time greatly over a fully two-dimensional method, the reason being that is 
much greater than the sound speed over the whole disk. Indeed, with this innovation we have 
been able to make hundreds of computer runs on a workstation testing various configurations. We 
have modified the method by Masset in several aspects for our codes and the full details of the 
numerical method are given in the appendix. 

2.4. Boundary conditions and Numerical Dissipations 

The boundary condition along the azimuthal direction is periodic. It proves to be very 
difficult in determining what is the best radial boundary condition. Since we are simulating a 
small part of the whole disk, ideally, we want to use the propagating sound wave conditions to 
minimize the possible reflection effects at the boundaries. Let 77 and T2 be the inner and outer 
disk boundaries, respectively. When "signals" are generated near tq and propagate with sound 
speed c,s both inwards and outwards, it takes roughly (r$ — r\)/c s and (r2 — to)/c s to reach the 
inner and outer boundaries, respectively. So, if Cg/v^ro) ~ 0.1, then after about 2 revolutions at 
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ro, signals would have reached the inner boundary. In five revolutions, signals would have reached 
a distance roughly 3 times ro outwards. Furthermore, the radial flows near both boundaries are 
likely to be subsonic, thus the boundary conditions could in principle affect the flow inside. We 
have tried several different choices with r\ = 0.2, 0.4, r2 = 2.0, 3.2, in various combinations. The 
dynamic behavior of the flow near ro (say [0.5 — 1.5]ro) is mostly independent of the size of 
the disk, as shown partly by the linear theory (Paper II) and the following simulation results. 
Using a large r2/ro ratio is relatively simple since the dynamics at large r is smooth and evolves 
slowly. Having a large ro/ri ratio, however, is obviously difficult. The strong Keplerian shear will 
continuously shorten the radial wavelength of radial propagating perturbations so that it becomes 
impossible to resolve them at late time. In most runs presented here, we have used r±/ro = 0.4 
and r2/Vo = 2.0, though we have also made many runs with ri/Vo = 0.2 and r2/Vo = 3.2 to ensure 
that similar results are obtained. 

After many tries, we have determined two types of radial boundary conditions that give 
reasonable results. One type of boundary condition is simply setting all the quantities in the 
ghost cells to be same as their initial values. This roughly mimics the condition that the two 
fluxes on the ghost cell boundaries are the same so that the mean quantities at the cell centers 
are not changed. Note that this condition still allows the material to flow off the grid since the 
flux is calculated at the cell interfaces. Another type of boundary condition for obtaining the 
ghost cell quantities is to fix and to extrapolate the variations between successive timesteps 
for the other three variables. The idea is to mimic the passage of "weak signals" assuming that 
the variations are typically small. We find that both types of boundary conditions work quite 
well on the outer boundary and there is negligible reflection from the outer boundary. At the 
inner boundary, however, we do not believe we have found a proper (or the most ideal) boundary 
condition, if there is one. The strong shear (i.e., short radial wavelength), the subsonic motion, 
and the fact that variations are incident with an angle into the inner boundary make it very 
difficult to eliminate reflections completely. Consequently, the mass flux through the boundary 
might not be quantitatively accurate, even though we have always observed mass flowing out of 
the grids in all our runs. 

Another important issue is the role of numerical dissipation. As discussed previously and 
further in the Appendix, we have employed a hybrid scheme that uses a first-order method (thus 
more dissipative) in regions with shocks and sharp discontinuities. The question is whether 
the transport we observed from the simulations is dominated by numerical dissipation or true 
physical effects. It is difficult to directly measure the amount of numerical dissipation. We have 
performed several tests. First, for a given initial equilibrium without the perturbations, we are 
able to evolve the disk to a time of > 50 revolutions at ro with the disk staying in the same 
initial equilibrium. For example, the maximum radial velocity (normalized by f^(ro)) found is less 
than 10~ 6 , compared to zero initially. Another set of tests we did was to look for convergence of 
azimuthally averaged quantities at different times using different numbers of grid points for the 
radial (nr) and azimuthal (np) directions, including nr = 100, 200, and 400, and np = 200, 400, 
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and 800. We typically find that one can use the 200 x 200 grid for a quick run out to 20 orbits at ro 
(~ 20 minutes on a Dec Alpha machine) with reasonably good results, and use the 400 x 400 grid 
for high resolution runs (~ 2.5 hrs out to 20 orbits). One further test was to compute the fraction 
of Lax-Friedrichs (LF) flux (first-order accurate) versus Lax-Wendroff (LW) flux (second-order 
accurate). This fraction is zero when shocks are not present (i.e., pure LW flux is used). For runs 
with strong shocks, this fraction can reach up to a few percent briefly and locally in the disk. 

We want to emphasize that the physics of the instability we will discuss later does not depend 
on the boundary condition, at least not critically. This is very different from some other global 
hydro instabilities, such as the Papaloizou & Pringle instability. The different treatments of the 
boundaries, however, do result in some minor quantitative differences, e.g., in estimating the 
transport efficiencies. Similarly, these estimates are affected by the numerical dissipation as well. 
Nevertheless, we are confident in the general physics presented here, such as how and where the 
angular momentum is transported and how the matter flows. But we are less confident in some of 
the exact numbers presented. Perhaps more sophisticated numerical schemes and better boundary 
conditions can improve this situation. 



2.5. Diagnostics 

The timescale of the simulations is referenced by the orbital period at ro (i.e., time t = 2ir 
stands for one revolution at ro). For most runs, we are able to run the simulations out to time 
t = 126 which is 20 orbits at ro- This translates into ~ 80 orbits at r\ = 0.4 and ~ 7 orbits 
at V2 = 2.0. This gives the system ample time for the nonlinear interactions to develop, since 
the linear growth stage usually takes < 8 orbits (see below). Depending on the amplitude of 
the maximum (v r ), this duration also spans a few local accretion timescales, allowing us to 
probe the accretion dynamics. As we will discuss later, 20 orbits are probably enough since the 
radiative cooling, which is not included here, is expected to play an important role after this many 
revolutions. 

Most variables are normalized by their values at ro, such as radius r/ro, angular velocity 
r2/f2fc(ro), growth rate 7/^fc(ro), and density £/£(ro). From these, the sound speed is normalized 
by the Keplerian velocity v<f,k(ro) and its value is typically ~ 0.1. The corresponding pressure is 
then ~ T,c 2 s /T ~ 6 x 10 3 for r = 5/3. So, if the velocity variations are close to 0.1, shocks are 
expected to occur. 

Besides displaying sequences of global 2D distributions of various physical quantities, the 
time-evolution of the disk can also be studied through various azimuthally averaged variables. 
One key quantity is the radial angular momentum transport due to the r — <j) component of the 
Reynolds stress (ESvrSv^), which is usually set equal to aP in the Shakura-Sunyaev formulation, 
where a is a dimensionless parameter characterizing the angular momentum transport efficiency. 
We will discuss this further in the presentation of our results. 
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3. Results: Comparison with the Linear Theory 

3.1. Confirmation of the linear theory 

The dependence of RWI on various initial equilibrium disks with different A and azimuthal 
mode number m has been given in linear theory in Paper II. In this subsection, we show that the 
linear theory is confirmed very nicely by our nonlinear simulations. For a given initial equilibrium 
with a specific A, we can use the eigenfunction of a specific unstable mode (i.e., 5T,, 5v r , Sv^, 8P 
from the linear theory) as the initial small amplitude perturbations. The linear theory predicts 
that the unstable mode should grow exponentially with a certain growth rate 7 and mode 
frequency u r . Furthermore, as shown in Paper II, the unstable modes are global, so that an 
exponential growth should occur throughout the whole disk. This behavior will certainly change 
when the nonlinear effects become dominant. 

The upper panel of Figure |I| shows the early evolution of radial velocities at 3 different radii, 
r/rg = 0.7, 1, and 1.3. The initial small perturbations are based on an eigenfunction with m = 5 
and A = 1.4 for an NGB equilibrium. The linear theory gives oj r =4.95 and 7 = 0.243. It is clear 
that the unstable mode is exponentially growing throughout the whole disk with the same growth 
rate and mode frequency. A rough estimate from this figure gives lo t ~ 4.92 and 7 ~ 0.24, which 
are essentially the same values as predicted by linear theory. In the lower panel of Figure [l] we 
present a run using an NGB equilibrium with A = 2.5 and m = 5. The predicted 7 is 0.61 and u> r 
is 4.9. The estimated values from the nonlinear simulation are again in perfect agreement. 

The nearly perfect confirmation of the linear theory also serves as a good test of our nonlinear 
code as it resolves the mode and captures its exponential growth. Similar results are obtained for 
other types of equilibria as well as different azimuthal mode number to, which we do not present 
here. 

3.2. Local Axisymmetric Instability vs. RWI 

One key result from the linear theory studies (Paper II) is that there exists a range of A where 
the disk is unstable to RWI but stable to the local axisymmetric instability. For A larger than 
certain critical value (which depends on the details of disk initial equilibrium), in addition to RWI, 
the disk has a small region where k 2 + N 2 is less than 0, where k and iV are the epicyclic frequency 
and the radial Brunt- Vaisala frequency due to the radial entropy variation, respectively. This 
makes the disk also susceptible to local axisymmetric instabilities according to the Solberg-Hoiland 
criterion (Paper II). In actual disks, the disk evolution will depend on which instability has a 
higher growth rate. We suspect that RWI will likely play an important role regardless. This is 
because the growth rate of RWI is usually already quite high (> 0.3Sl(ro)) for large A. Even if an 
axisymmetric instability grows first, since the quantity k 2 + N 2 is less than only in a very small 
region near rg (see Figure 5 in Paper II), the instability acts to stabilize this region, but might still 
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leave a finite A from which RWI can grow. Furthermore, since RWI is a global mode (compared 
to the axisymmetric instability which is local), its impact on the disk dynamics could be much 
larger. In some of our simulations that have a localized k 2 + N 2 < 0, it seems that RWI is always 
the dominant instability and in fact we have never detected any deviation from the exponential 
growth of RWI. We conclude that local axisymmetric instability is not important in our studies. 

4. Results: Nonlinear Stage 

Even though we are only dealing with 2D disk simulations, the nonlinear evolution of the 
flow is quite complicated. We will show in this section that large scale structures, such as vortices 
and in some cases, shocks, are produced. Besides the initial exponential growth due to the linear 
instability, even more localized pressure variations are produced both in the radial and azimuthal 
directions which feed back to the original instability. So, the disk evolution enters into a somewhat 
self-feeding state during which significant transport of mass and angular momentum is observed. 

4.1. Different Types of Runs 

We have performed a large number of runs using various initial equilibria. We choose to 
present 13 runs. Their properties are summarized in Table El. All 13 runs use nr xnp = 400 x 400. 
The first 12 runs have 0.4 < r/ro < 2.0 and are perturbed (from equilibrium) using their respective 
eigenfunctions with a specific azimuthal mode number m. All these simulations are run to time 
t = 126, i.e., 20 orbits at r = r^. The last run, T13, has 0.2 < r/ro < 2.0 with random initial 
perturbations. It is run to time t = 200 (i.e., ~ 32 turns). 

Runs T1-T4 represent 4 types of initial equilibrium with roughly the same linear growth rates 
(7 ~ 0.1) and similarly for runs T5-T8 (7 ~ 0.3). Runs T1-T4 all have small A so that n 2 + N 2 is 
everywhere positive (i.e., only RWI is present). Runs T5-T8, however, have a narrow region with 
k 2 + N 2 < 0, though RWI seems to be the only instability present. Figure § shows the evolution of 
pressure for runs T1-T4 and Figure || is for runs T5-T8. In order to make the pressure variations 
more clearly visible, we have actually plotted r 3//2 P for runs of NGB and NSJ types to take away 
the r -3//2 dependence in the background pressure (the second and fourth rows in Figures I and |). 
The left column shows the initial pressure distribution, the middle column is at the time when 
the linear instability just saturates (t = 3, and 7 orbits for lower and higher growth rate runs, 
respectively), and the last column is at time of t = 20 orbits. 

The initially axisymmetric pressure distribution has broken up and became nonaxisymmetric 
with distinct, organized regions, which turn out to be vortices. In addition, large scale spiral arms 
around these vortices have developed. It is clear that these nonaxisymmetric features, such as 
the hot spots in pressure, are quite persistent. Furthermore, runs T5-T8 are evolving at a much 
higher rate than that of runs T1-T4, especially in the nonlinear regime as well. This can be seen 
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Table 1: The initial setup of all 13 runs, which are categorized by homentropic or nonhomentropic 
Gaussian bump (HGB and NGB) and homentropic or nonhomentropic step jump (HSJ and 
NSJ). Different bump/jump amplitudes (A) and widths (Ar/Vo) are represented. All the runs 
except one (T13) have used eigenfunctions from the linear theory (Paper II) as the initial small 
amplitude perturbations, which are quantified by the azimuthal mode number m and growth rate 
7 (normalized by f2(ro)). Run T13 uses random initial perturbations. 



runs 


Type 


A 


Ar/ro 


m 


7 


Tl 


HGB 


1.12 


0.05 


3 


0.10 


T2 


NGB 


1.22 


0.05 


5 


0.11 


T3 


HSJ 


0.4 


0.05 


5 


0.11 


T4 


NSJ 


0.3 


0.05 


5 


0.11 


T5 


HGB 


1.35 


0.05 


3 


0.27 


T6 


NGB 


1.55 


0.05 


5 


0.32 


T7 


HSJ 


1.2 


0.05 


5 


0.29 


T8 


NSJ 


0.8 


0.05 


5 


0.30 


T9 


HGB 


1.6 


0.1 


3 


0.10 


T10 


HGB 


2.4 


0.1 


3 


0.32 


Til 


HGB 


1.17 


0.05 


3 


0.15 


T12 


HGB 


1.25 


0.05 


3 


0.20 


T13 


NGB 


1.51 


0.05 
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by comparing the middle and right columns in Figures ^ and |3[ Note that the pressure in the 
inner region (r/ro < 1) shows an overall increase as shown in Figure ^. This will be discussed in 
detail in later sections as we believe that this is a clear signature of overall accretion. 

Figures || and ^ also show that the dynamic behavior of the disk can be roughly divided into 
three regions: near r$ (say, r/ro ~ 0.8 — 1.2), inner (r/ro < 0.8) and outer (r/ro > 1-2) parts. 
Vortices and shocks form near tq, and they constantly generate waves that propagate towards 
both the inner and outer parts of the disk. (We have confirmed that these waves propagate at the 
sound speed.) These sound waves, being continuously sheared by the background flow, develop 
into spiral waves that might eventually lead to shocks. 

In all runs, the disk evolution can be roughly divided into three stages: an exponential growth 
of small amplitude perturbations, the formation of vortices, around which shocks are sometimes 
produced, and the global mass and angular momentum transport. The exponential growth phase 
has been discussed previously. We now discuss the rest of evolution in detail. 

4.2. Formation of vortices 

In this subsection we take a closer look at the formation of vortices. We specifically study 
two runs, Tl and T5, since they have relatively simple initial configurations, such as constant 
entropy. Figure ^| shows a global view of the whole disk with vortices. The pressure distribution 
is color-coded and the overlaying arrows map out the flow patterns around ro after subtracting 
v <b(fo,t = 0) (i.e., in the comoving frame that has a azimuthal velocity of K^o^t = 0)). The 
upper panel is taken from Tl at a time of 7 orbits and the lower panel is from T5 at 3 orbits. 
Both runs are initialized with the m = 3 unstable mode. Even though the Keplerian shear is still 
the dominant flow pattern, vortices are clearly formed in the flow. 

One dominant feature of these vortices is that the vortical motion is anticyclonic (the "spin" 
axis of the vortex is opposite to the disk rotation axis) and the vortex encloses a localized high 
pressure region. The nonuniform pressure distribution along the azimuthal direction (i.e., the 
—dP/d(j) term) is the main driving force in the formation of vortices. Such nonaxisymmetry 
grows out of RWI directly (see Figure 1 of Paper II). Taking the flow at r/ro = 1 from T5 as 
an example, in Figure ||, we have plotted the evolution of pressure P (upper panel), azimuthal 
velocity (middle panel) and radial velocity v r (lower panel) for times t = 0, 1, 2, 3 orbits, which 
are represented by the solid, dotted, dashed, and long-dashed lines, respectively. The flow is 
initially rotating with a single v<j> that balances the gravity plus pressure gradient. As the pressure 
gradient —dP/d(f> builds up due to RWI, the azimuthal velocity will decrease (increase) if 
the —dP/d(p force is negative (positive). A decreasing (increasing) causes the flow to move 
inwards (outwards) radially since gravity and rotation are no longer in balance. So, surrounding 
each localized high pressure region, an "anticyclone" is formed. In fact, this is probably the only 
vortical flow pattern that could survive in this nearly Keplerian shear flow. 
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The development of nonlinearity is clearly seen from these curves as well. At early times the 
azimuthal variations are still sinusoidal but become strongly concentrated by the time of 3 orbits 
(most clearly seen in the long-dashed curve for pressure, for example). Interestingly, even though 
the flow is nonlinear, these vortices remain azimuthally separated and are moving with the same 
speed around the disk (cf. Figures ^ and || at 20 orbits). Later we will discuss situations when 
vortex merge does occur. 



4.3. Vortex Radial Width and Shocks 

Figure ^ shows a closeup view of one vortex from those in Figure but now in the {r, </>} 
plane instead. The upper row is for Tl and the lower for T5 with pressure P on the left and 
entropy change AS = ln(P/S r ) — So on the right. Note that for both Tl and T5, the whole disk 
is set up with a single entropy So = ln(Po/^o)- 

As the flow is still predominantly Keplerian, there is a fundamental constraining effect on the 
radial width of a vortex. The Keplerian shear flow implies that the relative azimuthal flow speed 
(V(j)(r2) — v<f,(ri)) exceeds the local sound speed when 

\r2 ~ ri\/r > 2c s /v <t> (r ), (12) 

which is roughly ~ 0.2 in our case. In other words, imagine having a rod with a total length 
of 0.2ro and placing it radially with its mid-point at r$, the two ends of the rod could still 
communicate at the sound speed. This will not be true, however, for structures with larger radial 
width. In fact, this perhaps is the basis for the long-held belief that it is extremely difficult to 
maintain a long-lived vortex in accretion disks since it will be sheared away. 

The radial width of the vortices produced in runs Tl and T5 can be best estimated from 
Figure |^, where the streamlines are depicted for the same flow regions as in Figure ^. The vortex 



from Tl indeed stays within the limit imposed by equation (12) but the vortex from T5 is rather 
large. By "large" we mean that the vortex has a "core" region that has a radial extent of nearly 
QAtq. It also has extended "arms" (high pressure regions, cf. Figure @) that go out much farther in 
radial extent. The flow around this vortex is obviously much more complicated than that around 
the vortex from Tl. 

The critical difference lies in the fact that shocks are produced around the vortex in T5 but 
not in Tl. The vortex in T5 is surrounded by 4 "arms" which are labeled as A, B, C, and D in 
Figure fj] (see the corresponding locations in Figure ^). These arms mark the places where the 
pressure is high, so are the density and temperature (not shown here). Arms A and C are clearly 
shocks that have a pressure jump of nearly a factor of 2 and strong entropy production (lower right 
plot of Figure |6|). Note that arms A and C start at the radial location r ~ 0.9, l.lrp, respectively, 
agreeing precisely with that expected from equation (|i"2|). Arms B and D are, however, not shocks. 
Instead, they are the location where the rarefaction wave from shocks at arms A and C meets with 
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the background flow. This is supported by the fact that no entropy variation is seen even though 
the pressure variation is close to a factor of 1.5. In addition, the streamlines are smooth near arms 
B and D (unlike A and D), the apparent flow direction change at arms B and D are actually due 
to an imbalance of gravity, rotation and pressure forces, not from a shock. The flow is expanding 
very strongly in the region between arms B and C, as well as between A and D, where the flow 
directions are strongly altered and the radial velocity reaches its maximum for both the infalling 
and outward motion (indicated by the length of the arrows in Figure ||). We believe that these 
flow structures, especially the shocks, are pivotal in the formation and "protection" of the vortices 
against the background shear flow. They have enclosed a region within which the flow is subsonic 
and the streamlines are closed. 

In contrast, no shocks are produced around the vortex in Tl, as evidenced by the extremely 
small change in entropy shown in the upper right plot of Figure ||. (Note the different scales of 
AS in the two entropy plots.) These changes at ~ 10~ 4 level are most likely due to the numerics 
alone. In other words, the whole disk remains homentropic to a high degree. 

4.4. Dependence on the initial width At/tq 

We have also investigated the vortex radial size dependence on the initial pressure bump 
width Ar/rQ. This is done in runs T9 and T10, where Ar/r^ = 0.1. These two runs are designed 
to have the same linear growth rates as runs Tl and T5, respectively. A very similar evolution is 
observed in these two runs, and Figure |8] shows the flow velocities around a single vortex together 
with its pressure distributions in the {r, 0} plane for runs T9 (left) and T10 (right). They are 
taken at times of 7, 4 orbits, respectively. The lower panels are shown at a time of 20 orbits. 

The vortex in T9 has a small radial width (< 0.2ro)), the same as that in Tl. Consequently, 
no shocks are observed either. For the vortex in T10, shocks are again observed at very similar 
locations as those in T5. This can be seen by comparing the upper right plot of Figure |8| with 
the lower left plot of Figure ||, the two shock structures (arms A and C) are quite similar, with 
the same starting radial locations at ~ 0.9 and l.lro, respectively. The compressions due to 
rarefaction waves, however, are not as pronounced in T10 as those in T5. Thus the radial width of 
a vortex does not depend on the initial pressure bump width. The characteristics of the Keplerian 
flow is the dominant factor in determining the vortex width, with or without shocks. 

4.5. Radial Drift of Vortices 

The right panel of Figure || also reveals an interesting phenomenon: there is a slight but 
noticeable inward radial drift of the vortices between the times shown (4 and 20 orbits). At 20 
orbits, this drift is only visible in run T10 but not in T9, presumably because T10 is evolving much 
faster. The amount of radial drift appears to be small, but it actually implies a high accretion 



-15- 



rate. Using the usual scaling relation for the radial accretion velocity, v r ~ ac s (H/r), where H is 
the disk vertical scale- height, we get 

v r r Ar drift 2 
a Z „ . (13) 

ity c s H 2nr N c s 

where N is the number of orbits at rg, and v^/cg ~ 10. Reading from the right column of Figure 
||), we get Ardrift ~ 0.05ro and N = 20, implying that a ~ 0.04. This simple estimate turns out to 
be quite consistent with more detailed analysis presented later. 



4.6. Dependence on the initial perturbations 

In physical systems, the initial perturbations are unlikely to be a single eigenfunction given 
by the linear theory, though one can obviously decompose the variations into various eigenmodes. 
In run T13, we use a random small amplitude initial perturbation (this ensures nonaxisymmetry 
by default). An initial exponential growth is again observed (not shown here). In Figure ^, 
we show 12 snapshots of the disk in color-coded radial velocity, with the first 11 frames at 
t = 0,2,4,6, ...,20, and the last frame at t = 32 orbits. Vortices are clearly produced, just as in 
all the other runs we have presented. There is a clear trend that vortices merge with each other, 
going from ~ 5 — 6 vortices initially to only one dominating vortex at ~ 16 orbits. The fact that 
the m = 5, 6 modes grow first is expected from the linear theory analysis as they have the highest 
linear growth rates (cf. Figure 10 of Paper II). Note that even though the vortices are nearly 
corotating with the background flow, there is nevertheless a difference in the phase velocity for 
different m modes. Eventually a faster moving vortex will catch up with a slower one and the two 
vortices will merge. In the end, there is only one vortex left in the system since, given enough 
time, any slight difference in the phase velocity will lead to an interaction between two vortices. 
We emphasize that such strong interactions are due to the fact that vortices are excited/produced 
at nearly the same radius and the radial drift of these vortices are very slow. In a real system 
where multiple "bumps" might be present at different radii, multiple vortices could be present. 



4.7. Mass and angular momentum transport 

We now address the critical question of the mass and angular momentum transport in these 
disks. Neglecting dissipation (e.g., shocks) for a moment, angular momentum is conserved (except 
for the loss due to the flow through the boundaries), but may be redistributed as the disk evolves. 
Following the treatment in Balbus & Hawley (1998), we can separate the velocities into a mean 
component and a "varying" component, v^t) = v^t) + Sv^ and v r (t) = + Sv r (t), where v^it) 
is obtained by averaging v^r, <j>, t) over <j> at a particular radius r and time t. Consequently, the 
radial flux of angular momentum is decomposed into two parts 



r 2 [v^(Ev r ) + (EvrSv^)] , 



(14) 
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where (...) indicates averaging over f d(j)/2ir. The first term indicates the direct radial flow of 
matter and is proportional to the mass accretion rate, 27rR(T,v r ). The second term represents the 
radial angular momentum transport through the correlations of velocity component variations. 
Traditionally this has been thought of as the turbulent Reynolds stress (T r ^) = (TiV r 8v^)^ whose 
origin has been the subject of intensive research for decades. Furthermore, as emphasized in 
Balbus & Hawley (1998), what is more important is the positive correlation between v r and Sv^ 
instead of the mere presence or amplitude of these variations. Even though we do not regard the 
flow we are studying as turbulent, the same requirement, i.e., the positive correlation between v r 
and <5iu, still holds the key to an outward transport of angular momentum. 

The most important result of the vortices generated by RWI is that the flow pattern around 
these vortices is perhaps an ideal configuration for an outward angular momentum transport 
process. This is due to the fact that the azimuthal pressure gradient causes variations in and 
consequently leads to the generation of v r via radial force balance, i.e., a decrease (increase) in 
Vff, leading to a negative (positive) v r (see detailed discussion in §|4^). Such a correlation directly 
ensures the radial angular momentum flux via transport (cf. equation ([lj)) is positive, i.e., an 
outward transport of angular momentum. 

To quantify the crucial role of vortices in angular momentum transport, we define a 2D 
version of a modified a coefficient, 

_ Sjj v ri j [v^ij -Ty,(r,t)] 

where the indices {ij} stand for {r<p}. Here, we use v^r, t) as the "mean" background flow (though 
it is a debatable choice), but it nevertheless gives a good indication as to which regions/structures 



are contributing most importantly to the angular momentum transport. Figure 10 shows the 
strength and distribution of a r( p around a vortex in the {r, <p} plane for run T5 at time of 3 (left 
panel) and 20 orbits (right panel). Similarly, Figure |ll| shows the distribution of ct r< ^ at 7 (left 
panel) and 20 orbits (right panel), respectively, for run Tl. 

Comparing Figures |(| and ffoj, one can see that strong outward angular momentum transport 
occurs in the expansion regions behind shocks A and C, with localized a r( j, exceeding 0.1. Similar 
structures are observed in results from Tl even though the amplitude is at a much reduced level 
and shocks are not present. Even after 20 orbits, the main features and their strength (notice the 
scaling) remain amazingly steady and strong. These vortices are rather remarkable in this regard. 

To further quantify the global transport efficiency, we can take an azimuthal average of 
equation (|l~5|), which is equivalent to 

(a) = (T r(j> /P) , (16) 

where {T r( p) = (Edvrdv,),) and 



(ESvrSv,/,) = (Ev^) - (T,v r )(£)/r 



(17) 
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where (I) = is the average specific angular momentum (Hawley 2000). The lower panels 

of Figures [l0| and |ll] give (a) at the same times as their corresponding upper panels. The transport 
efficiency of T5 is much higher than that of Tl, by a factor of at least 30 (though there is only a 
factor of 3 difference in linear growth rates). So, it is not surprising that the evolution of T5 is 
much faster than that of Tl, as observed previously. 

To illustrate the dynamics of transport efficiency, we present the evolution of (a) for runs 
T1-T8 in Figure |l2|. The strength of transport usually reaches a peak when the vortices first form 
(the dotted lines), but settles down to maintain a steady level, and so there is relatively little 
difference between 10 and 20 orbits. Even though the transport of angular momentum in the disk 
has both outward and inward components (as indicated by the positive and negative values of a r ^ 
in Figures and 0), on average the angular momentum is transported outward through each 



radial "ring", as indicated by the predominantly positive a given in Figure 12 for all the runs. 



This is extremely encouraging and perhaps the most important result of this study. 

Another important point is that the transport can be roughly divided into two different 
regions, that associated with vortices and that associated with the trailing spiral waves that are 
present in both the inner and outer parts of the disk. The physics behind the outward angular 
momentum transport by trailing spiral waves is actually quite similar to what we have discussed 
previously for the transport by vortices, since the azimuthal pressure gradient is fundamentally 
responsible for causing the positively correlated velocity variation components. In some sense, a 
vortex is just a much more pronounced nonlinear manifestation of such transport processes. In 
the case that the vortex is radially large and strong (i.e., run T5), the transport strength around 
the vortex is much larger (a ~ 10~ 2 ) than that of the spiral wave region (a ~ 10~ 3 ). In the case 
that the vortex is weak (i.e., run Tl), they become comparable («~3x 10~ 4 ). The increase of a 
at smaller radii is related to the stronger shear, though it is difficult to accurately estimate it. 

Note that the spiral waves are produced by continuously shearing the radially propagating 
sound waves that are generated in the vortex region. This is one of the important features of the 
linear theory (Paper II) in which waves are allowed to "tunnel through" the trapping region where 
the vortex is produced. The consequence of this connection is that transport will occur not only 
near the position of the initial bump/jump, but also throughout the disk as a whole, thus giving 
rise to a much larger, global impact. 

We have also investigated other physical quantities using azimuthal averages. These quantities 
include the pressure (P), mass surface density (£), mass flux (accretion rate) F p = J dcprYlvr, and 



angular momentum flux Fj = J d^rJ^rv^. They are shown in Figures 13 and 14, for runs Tl, T6 
and T8, respectively. 

There are several generic features that appear in all these runs: 

(1) As expected, the instability always tries to remove the bump/jump. The lower linear 
growth rate runs evolve more slowly than those with higher linear growth rate runs. The nonlinear 
saturation levels are also different, as evidenced by the magnitude of mass accretion rate F p , for 
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example, in Figures 13 and 14 



(2) As a clear confirmation of the efficient accretion that is going on in the disk flow, 
the pressure in most parts of the disk is increasing with a variable amount (except the initial 
bump /jump, of course). This is especially true for the inner part of the disk. We believe that this 
increase is from the release of gravitational energy due to the "global" accretion caused by both 
vortices and spiral waves. This effect seems inevitable since we did not include any cooling effect 
in our equation of state. 

(3) The global nonlinear evolution brings an additional lack of axisymmetry and radial 
variations in the disk flow. This is manifested in the average pressure distributions where large 
radial as well as azimuthal (cf. Figures [2] and |||) gradients are produced. These strong gradients 
will be susceptible to the same Rossby wave instability we are studying. It is then not very 
surprising that the system can sustain itself for a long time, consistent with our results. 



5. Discussion 



One advantage of our nonlinear simulations is that we are guided by a robust linear instability 
that has been investigated previously. The precise confirmation of the linear theory not only 
validates the presence of this instability, but also provides a check for our nonlinear codes. 
Consequently, some of the usual concerns with numerics are not as important. At the fully 
nonlinear stage when shocks are present, it is, however, difficult to capture all the dissipation 
perfectly. So we are less confident in some of the exact numbers presented, but we believe that the 
large scale structures of this instability have been captured correctly. Furthermore, by following 
the instability evolution through the linear growth stage, we have gained more confidence on 
the physical mechanism of the instability and have singled out the key physical processes in the 
nonlinear regime, such as the formation of vortices and shocks. There are still a number of physical 
issues that deserve further discussion. 



5.1. Setting Up the Initial Equilibrium 

In realistic astrophysical situations, the initial conditions will certainly be system dependent. 
The idealized bump/jump along with the background disk described by the present studies might 
arise in the close binary systems where matter tends to be stored at the large radii first; Or it 
could be the radiation heating from the central star that causes a localized region of the disk to 
be hotter; Or it could be the edge between a proto-planet and its surrounding disk material. By 
investigating the role of RWI in many different initial configurations, we show that RWI is very 
robust. Quite generally, the disk is potentially unstable due to RWI whenever there are "bumps, 
edges, clumps" present in the disk. 
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5.2. Dependence on the amplitude A and Physics of Saturation 

Since there is a clear difference between runs with low and high growth rates, we now compare 
4 runs, Tl, Til, T12 and T5, all of which use a Gaussian bump but with an increasing amplitude 
A. Their linear growth rates range from 0.1 to 0.27. Instead of showing global distributions of 
various quantities, we opt for a single quantity, the maximum radial velocity \v r (ro, <fi, t)\ at ro, 
to describe the development of the instability. This is plotted in Figure [l^ for different runs as 
a function of time. Along with these 4 runs, we also plot the results from runs T6, T7, and T8. 
While it seems that all runs with high growth rates eventually saturate at roughly the same 
max(|tv|), the saturation level of max([tv|) for lower growth rate runs clearly depends on the linear 
growth rate. 

Since the value of max(|u r |) can be regarded as a rough measure of the level of nonlinearity 
and the associated transport efficiency, one question that naturally arises is "what is the physics 
causing the instability growth to saturate?" There are several possibilities. One is that when 
the linear growth is slow, saturation is achieved by the removal of the bump /jump, i.e., the 
driving of the instability. We believe this is what happens when the bump/jump is small, such 
as runs T1-T4. The evolution shown in Figure [l3| supports this interpretation (see also Figure 
H). As shown before, anti-cyclonic vortices are formed surrounding regions with high pressure and 
density, but there are no shocks in the flow and the flow entropy is well conserved (cf. Figure ||). 

A relevant study on this issue can be found in Laughlin et al. (1997, 1998) with a somewhat 
different setting. From their excellent nonlinear mode coupling analyses in self-gravitating gaseous 
disks, they concluded that the growth of the dominant unstable mode can modify the background 
disk profile so as to prevent its further growth, causing saturation without relying on dissipation. 

When the bump/jump gradually increases, the instability growth becomes fast enough that 
the saturation is achieved both by the removal of the bump/jump and the formation of shocks. 
Anti-cyclonic vortices are again formed with their radial width being roughly 4 times the thickness 
of the disk. The larger radial size leads to the formation of shocks, which in turn limit the growth 
of radial velocities, causing the instability to saturate. 

5.3. Energy Conversion and Dissipation 

One of the key physical links in accretion disk physics is to understand how gravitational 
energy is converted into internal energy, part of which can then be radiated away. The increase 
of internal energy, however, can occur either adiabatically (i.e., entropy is conserved) or with 
dissipation (i.e., heat generation with increasing entropy). In the classical, geometrically thin and 
optically thick a disk model (similar to the physical condition we are considering here), the disk 
is assumed to be axisymmetric, Keplerian and quasi steady. The angular momentum transport is 
via an assumed anomalous viscous stress that is related to the Reynolds stress (Su r (5^}. These 
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conditions lead to a relation where, at large radii, local heat production (via dissipation) can be 
a factor of 3 times the available gravitational energy release, a direct result of viscous transport 
(cf. the textbook by Frank, King & Raine 1985). The viscous heating rate per unit volume per 
gram is dQ ~ v(rdQ/dr) 2 (Lynden-Bell & Pringle 1974). Estimating v ~ ac s H, where H is the 
thickness of the disk, we would expect an entropy increase dS by relating TdS = dQ, 

dS (9/4) {fit) a , (18) 

i.e., the entropy increases linearly with time as transport continues. 

Our results in run Tl, however, seem to indicate a different route in the energy conversion 
process. The evolution in Tl satisfies the adiabatic condition to a high degree, as evidenced by 
the near constancy of entropy of the whole disk, which maintains its initial value (cf. Figure |6|). 
(In fact, we have written a different version of the code by requiring that the entropy of the flow 
be conserved, as contrasted with the code we presented which uses the total energy conservation. 
Both codes give very similar results.) This disk, however, is transporting angular momentum 



outward with an equivalent a of 10 -4 (conservatively), as shown by Figure 11. Using equation 
fli~8D, this level of angular momentum transport would imply an entropy increase far greater than 
what we have obtained, which is less than dS ~ 2 x 10~ 4 as shown in upper right plot of Figure^. 
This special case proves an important physical point, that angular momentum can be transported 
outward in a disk without dissipation. The released gravitational energy goes entirely to PdV 
work, which is done adiabatically. Again the studies by Laughlin et al. (1997, 1998) are relevant 
here. In their case, the nonlinear mode interactions create a nonsteady perturbing potential that 
continuously drives the disk evolution, giving angular momentum and mass transport without 
dissipation. 

Whether the real astrophysical disk operates via a "maximal heat/entropy production" route 
or an adiabatic route is unclear, and observational constraints have been scarce. Different energy 
conversion processes, however, do predict different amounts of energy that can be radiated and 
consequently give different radiation spectra. It is needless to say that the global disk structure 
and evolution obtained from our study differ fundamentally from the classical a disk model, 
especially in how energy is converted and how heat is generated. For example, one obvious 
difference is that the transport efficiency has a radial dependence, which implies that the disk is 
always in a dynamic state and can be more properly described as "surges" of matter accretion. 
How to meaningfully construct global accretion disk models under these conditions and relate 
them to observations might be very interesting. 



5.4. Late Time Evolution and Effects of Radiative Cooling 

We are able to run most simulations out to 40 — 50 revolutions (at 7"o) and find that disks are 
continually evolving, though we have only presented the results up to 20 orbits. In fact, in some 
cases the disk has evolved so much that we believe that we should not run those simulations much 
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longer than 20 revolutions. This is because additional physical effects that are not included in the 
present study could become too important to ignore. One such effect is the radiative cooling. For 



example, as shown in Figure 14, the average inner disk pressure has gone up by a large factor (i.e., 
x2), while the density shows only relatively small variations. This implies a large change in disk 
temperature, which could mean a large change in radiative loss as well. 

As emphasized in Paper II, the disk has to be relatively hot (i.e., c s jv$ > 0.05) in order to 
have a "healthy" growth rate for RWI. The increase in disk pressure helps to sustain the dynamic 
evolution. On the other hand, ignoring radiation means that the cooling time of disk should be 
relatively long compared to the disk rotation periods. This requirement, as discussed in Paper II, 
implies a minimum column density of the disk matter so that heat can be trapped inside the disk 
for several revolutions. Since our simulations are 2D, we could not directly model the effects of 
radiative cooling with respect to the vertical transport. One way to circumvent this difficulty is 
to add an ad hoc local cooling function that removes the internal energy at a specified rate (see 
Flozyczka & Spruit 1993). Similarly, we have tried to add a loss term in the energy equation as 
— e/r c , where e = P/iT — 1) is the internal energy and r c is the characteristic cooling timescale. 
The parameter r c is likely to be a complicated function of density, temperature (hence radius of 
the disk) and radiative transport processes. As a simple approximation, we have tried to relate 
r c to the local Keplerian rotation period by a single constant. Indeed we find the trend that the 
shorter the cooling time, the weaker the nonlinear effects of the instability. Intuitively, if the 
cooling time is shorter than one rotation period, then changes caused by the vortex motion (i.e., 
compressions and expansions occurring during a "thermal" cycle) will likely be damped out very 
quickly, thus strongly limiting the efficiency of any transport processes. 



5.5. 2D versus 3D 

Another important aspect is the 3D nature of the disk flow. The 2D approximation is 
expected to break down in several ways. The effects of strong shear are clearly visible in all the 
runs (cf. the inner region of Figures ^ and |3|) where the spiral waves become wound tighter and 
tighter with increasing orbital velocity. Such a short radial wavelength situation probably violates 
the 2D assumption. The radial propagation of the sound waves itself can be weakened by the 
"refraction" effect discussed in Lin et al. (1990), thus limiting its impact radially. Furthermore, 
as the disk expands/contracts both vertically and radially as pressure varies, dissipation by the 
irreversible processes becomes inevitable (e.g., expansion into a near vacuum). This will probably 
prove to be the most important 3D effect, though detailed 3D simulations are needed to address 
this problem quantitatively. 
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5.6. Entropy gradient is not necessary 

It is worthwhile emphasizing that RWI grows under a wide variety of physical conditions and 
that it is not necessary to have an entropy variation in the disk, at least for 2D disks (Paper II). 
This is supported by the above results (run Tl, for example). What drives the mode unstable is 
the steep pressure gradient which gives rise to the "trapping" physics that allows the mode to be 
amplified. The amplification of the unstable mode is the result of repeatedly passing through the 
corotation radius that is residing in a "trapping" region (Paper II). 

Entropy variations in the disk could, however, introduce more features since the additional 
potential vorticity can be driven thermodynamically by the VP x VS term. In actual disks, 
entropy variation might be inevitable as argued by Klahr & Bodenheimer (2000). Further studies 
are needed in order to better quantify the effects of entropy variations. 

6. Conclusions 

We have studied the global nonlinear evolution of the Rossby wave instability in a nearly 
Keplerian flow, following our previous linear theory analysis (Lovelace et al. 1999; Li et al. 2000). 
During the linear growth stage of the instability, our nonlinear simulations agree extremely well 
with the linear theory results (e.g., the growth rate and mode frequency). 

In the nonlinear stage, we have shown that vortices naturally form, enclosing a high 
pressure and density region. These vortices are (nearly) corotating with the background flow 
but are counter-spinning (i.e., anti-cyclones). We have elucidated the physical mechanism for 
the production of these vortices, namely through the azimuthal pressure gradients, and shown 
that they are long-lived structures within disks. These vortices are shown to be extremely 
important for transporting angular momentum outward and for causing global accretion. In fact, 
by analyzing the flows around each vortex, we have shown why they are the ideal "units" for 
outward angular momentum transport, namely by giving rise to "positively" correlated radial and 
azimuthal velocity variations, i.e., (YiV r 5v<f,) > 0. Shocks are formed when the instability is strong 
and these shocks limit the radial extent of a vortex to be less than 4 times the thickness of the 
disk. Furthermore, trailing spiral waves are produced both interior and exterior to the vortex 
region. These trailing spiral waves are produced by shearing the radially propagating sound waves 
generated by vortices, and they serve as an additional means of transporting angular momentum 
outwards. 

We find that the angular momentum transport efficiency is not a constant throughout 
the whole disk; it has both radial and azimuthal dependences, and evolves continuously. 
Consequently, we recognize the need to construct global models of accretion disks that reflect such 
nonaxisymmetry and dynamics. 

Several important physical issues are, however, not addressed in this paper. Three- 
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dimensional modeling is especially needed to address the issue of radiative cooling and how much 
heat/dissipation is produced locally. In addition, combining this hydrodynamic instability with 
those associated with magnetic fields will be very interesting as well. 
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A. Numerical details 

Our numerical method differs from that used by Masset (2000) in several ways: 

(1) We do a Galilean transformation of the split angular equations, transforming to a 
coordinate system rotating with a constant velocity at each radius. The velocity is chosen to be 
as close as possible to the mean angular velocity such that the transformation back to the fixed 
system involves only a shift of indices. 

(2) The stability condition is the standard Courant, Friedrichs, and Lewy (CFL) condition 
computed from the radial velocity and the sound speed, with limit 1. The time interval determined 
in this way is still too large for the angular integration, in spite of the reduction in angular velocity 
obtained from the co-moving system, but we use partial time-stepping in angle, satisfying the 
angular CFL condition at each partial step. 

(3) All the sources are included in the radial sweep; the sources are not done in a split step. 

(4) A simple hybrid scheme is used for both the angular and radial integrations. The idea is 
to use a weighted combination of the second-order Richtmyer two-step version of the Lax-Wendroff 
(LW) scheme and the two-step first-order Lax-Friedrichs (LF) scheme, with weights chosen to favor 
LF in regions possibly containing shocks, and to favor LW in smooth regions, while maintaining 



the second-order accuracy there. This is an idea first proposed in Harten fc Zwas (1972) . There 



are many variations of this method that are described very well with references to the original 



literature in Laney (1998). The particular weights we have chosen seem to provide a robust 
scheme. 
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In dimensional splitting we consider separately the radial equations with the sources; 



and the angular equations 



-^w + —f(r,w) + S(r,w)=0, (Al) 



m w + ^ w) = - (A2) 



A single cycle consists of first a determination of the time interval At, then a sequence spanning 
two time steps of the form radial-angular-angular-radial. More formally, we would have calls to 
routines of the form 

CALL radial(w,w') 

where the input w is the vector of dependent variables at time t n , and w' is the output. Then 
using w' as input, 

CALL angular (w',w) 
where w is now the result of one time step. 



A.l. The radial integration 

Equal radial intervals Ar define equally spaced radii rj, i = 0, N + 1, presumed at the center 
of radial cells. The cell vertices are at r i+ i . The inner radial boundary is at ri , the outer at r N+ i. 

In describing the radial integration we suppress the dependence on angle: the difference 
equations being described must also be applied at each angle. The equations are written in flux 
form 

w[ = Wl - \{F l+l - F t _i) - AtS*, (A3) 

A = Ai/Ar, (A4) 

with the flux F and the source S* defined below. 
The flux is a hybrid: 

F = aF LF + (1 - a)F LW , (A5) 



< a < 1. 
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F LF is a Lax-Friedrichs flux: 



i£i = -^K+i - ^) + + /K)) + (A6) 



(A7) 



F w is the Lax-Wendroff flux 



«+7 



(A8) 



The source S* is a vector with four components, S rn *,m = 1,4, but only two are nonzero, 
S 1 * = S 3 * = 0: 

4* / r*4* i o4* \ 



with 



and, 



Sf* = -(S 4 *! + s 



5 



■4* 



(A9) 



(A10) 



-,2* 



+r i+ i(P(w, 



l/r)* + i + (£((^) 2 
)- j P(<l))/dr. 



l/r)J 



(All) 



Turning to the predictor step equation( [A7|) , we found it necessary to make an adjustment. The 
greatest stumbling block to a long run, say ten turns of the disk, is the appearance of a negative 
pressure somewhere in the disk. The final set of difference equations we have used does not have 
this defect without setting any artificial lower bound on the pressure, at least in all the many runs 
we have done. The problem is that apparently minor changes in the scheme can cause negative 
pressures. Thus, we found that using the angular momentum conserving form in the predictor step 
did not work. Of course, angular momentum is still conserved, so there should be no objection on 
that ground. 



The weight a is 



o 



min( 



(V r i-l 



Ci 



(A12) 



where q is the local sound speed. 
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A.2. The angular integration 

The angular integration finite-difference equations have the same hybrid form as those for the 



radial integration, except that there are no sources, and replaces v r in ( A12 ). Since the radius 
is constant on any one angular sweep, all extraneous radial factors can be canceled. The hybrid 
weights are computed using angular velocity differences rather than radial velocity differences. 
However, as indicated in section |2.3| , a local co-moving frame and partial time-stepping are used. 

Equal angular intervals A(f> define equally spaced angles (f>j,j = 1,M, presumed at the center 
of angular cells. The cell vertices are at 6- , 1 , with 6\ = 0, d> M ,i = 2tt. The periodicity is enforced 

JT 2 2 ' 2 

by means of the indexing, that is, the index j is always understood as j + M — int((j — 1+M)/M)M. 

In describing this method, we suppress the radial dependence. The local co-moving frame is 
determined by first defining the average angular velocity, following Masset. 

^ = ^£f(Vb'- (A13) 



The actual velocity shift used is given by 



Vcj, At 
r A(f> 



(A14) 



Z = int(r + i) (A15) 

v = Ir^t (A16) 

The purpose of this shift is to reduce the effect of the large angular velocity on the time 
step for the angular integration. If the angular velocity were independent of j, the use of 
would accomplish this at the expense of having to do an interpolation later, since the shift of the 
coordinate system in one time step would not be an integer number of cells. By using vq it will 
be an integer shift Then, at the beginning we define new conserved variables by replacing tu with 
v '<t> = v <t> ~ u o- This entails a replacement of Hv^ by and of HE by jAj- + ^£((t^) 2 + v^.). At 
the end of the angular sweep we will have computed new variables, Sj, (v r )j, {v^)j, and Pj, but 
these are in the moving coordinate system. To get back to the fixed system, first replace (v^j by 
{ v <j>)j + u 0) then shift the density, velocities and pressure from cell j to cell j + I, then redefine the 
conserved variables. 

For the actual integration we cover one time interval At with an integer number of steps with 
interval 

At' = — (A17) 

n 
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where n is determined from the angular stability condition. Thus, if 



E = max( 



v '<p + c \ K~ c 



(A18) 



r 



r 



then 



n = int( 




) + l 



(A19) 



We have found that the average n is about 1.2. 
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Fig. 1. — The exponential growth of the radial velocity during the linear growth stage of the Rossby 
wave instability. Two runs, both of NGB type, are shown with different bump amplitude A. The 
dotted curves show the magnitude of radial velocities at three fixed locations in the observer's frame 
at r = 1,0.7, 1.3 and = 0, from top to bottom, respectively. For a given A and m (= 5 for both 
panels), the linear theory of RWI predicts a specific mode frequency uj r and growth rate 7. All 
the curves are from our simulations and from each curve, one can get u r from the oscillations (the 
dotted line) and 7 from the slope of its "envelop" (the solid line). Both quantities show excellent 
agreement with the predictions of the linear theory. 
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Fig. 2. — The evolution of the pressure for runs T1-T4 (lower growth rate runs). Each row consists 
of snapshots of the whole disk at three different times of each run (T1-T4 from top to bottom). 
From left to right, t = 0,7, 20 orbits, respectively. The color code is for pressure, which is in units 
of 10~ 3 . Note that the scale is different for each run. Each run is initialized using small amplitude 
perturbations based on the eigenfunction of its linear instability with a specific azimuthal mode 
number m, which is m = 3, 5, 5, 5 from top to bottom. The pressure of T2 and T4 (second and 
fourth rows) has been multiplied by r 3 / 2 in order to make the pressure variations more easily visible. 
Isolated hot spots (high pressure) are clearly visible, and they are the centers of large anticyclonic 
vortices. Large scale spiral arms are produced in connection with these vortices as well. 



Fig. 3. — Similar to Figure || except using results from runs T5-T8 (high growth rate runs). The 
middle column is now at a time t = 3 orbits and the right column is again at 20 orbits. The 
amplitude of the pressure variations is larger than that seen in Figure |2|. The pressure has clearly 
increased in the inner part of the disk. Vortices and large scale spirals are produced as well. 



Fig. 4. — Vortices in a disk. Pressure is color-coded (in units of 10~ 3 ). Arrows indicate the flow 
pattern near tq in a comoving frame of vAtq). Vortices are anticyclonic, enclosing high pressure 
regions. Large-scale spirals are produced as well, in connection with the vortices. The upper panel 
is a snapshot from Tl at t = 7 orbits and the lower panel is from T5 at t = 3 orbits. Both runs 
use an m = 3 unstable mode, which is why there are three (nearly) corotating vortices. 



Fig. 5. — The production process of vortices. Shown is how the pressure (upper panel), the 
azimuthal velocity (middle panel) and the radial velocity (lower panel) vary along the azimuthal 
((f)) direction as the instability develops. The results are from run T5 and are taken at r = r$. The 
solid, dotted, dashed and long-dashed curves in each panel are at t = 0, 1, 2, 3 orbits, respectively. 
Note the tt/6 shift between the peaks in pressure and peaks in v<f>. Such correlation is derived 
from the fact that the azimuthal pressure gradient —dP/d<f> is responsible for the variations in 
Va- Consequently, the imbalance between gravity, rotation and radial pressure gradient in the 
radial direction introduces radial motion. This explains why the largest positive (negative) radial 
velocities occur when is the largest (smallest), i.e., super- (sub-) Keplerian. 
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Fig. 6. — A closeup view of one vortex in the {r, (f\ plane. The upper two plots show one vortex 
from run Tl at t = 7 orbits. The flow near vq is indicated by the arrows, together with the 
color-coded pressure (left plot) and the color-coded entropy change AS (right plot). Similarly, the 
lower two plots show one vortex from run T5 at t = 3 orbits with the pressure (left) and entropy 
change (right) distributions. The vortex from Tl is relatively weak with small radial motions. No 
shocks are present, as shown by the exceedingly small variations in entropy AS. The vortex from 
T5, on the other hand, is very strong with significant radial motion. There are 4 pressure "arms" 
(locations with high pressure) surrounding this vortex, as indicated by the labels A, B, C, and D. 
Shocks are formed at arms A and C, as shown by the large increase in entropy. Arms B and D 
are probably not shocks, instead, they are produced by the rarefaction waves from the shocks at 
A and C. Consequently, entropy does not change at B and D. The high pressure band and the 
high entropy band in the lower left region of the two T5 plots (lower panel) are from the arm A of 
another vortex (not shown). 

Fig. 7. — The streamlines around the same vortices from Tl (left) and T5 (right) as those shown 
in Figure |6|. The vortex from Tl has a small radial width (< 0.2ro) but the vortex from T5 has a 
large radial width (~ 0.4ro). The sharp changes in streamline direction at arms A and C further 
prove that they are indeed shocks, whereas arms B and D are not. 

Fig. 8. — A closeup view of vortices from run T9 and T10, both of which have a wider initial 
bump width Ar/Vo = 0.1. The pressure is color-coded in all plots and the comoving flow velocities 
are overlayed as arrows. The two left panels are from T9 at t = 7 (upper) and 20 (lower) orbits, 
respectively. The two right panels are from T10 at t = 3 (upper) and 20 (lower) orbits, respectively. 
Again, shocks are formed around the vortex from T10, limiting its radial width to be less than 
~ 0.4ro. There is a noticeable inward radial drift of the vortex from T10. A high angular momentum 
transport efficiency is implied from such a drift. 

Fig. 9. — Merge of vortices. Shown are 12 snapshots from run T13, where random initial 
perturbations also lead to large-scale vortex formation. From top left to bottom right, the frames 
are from time t = 0,2, 4, 6, 20 and t = 32 orbits, respectively. The radial velocity is color-coded 
and the same color scale is used for all frames. The amplitude of the radial velocity grows from 
small values (~ 0) to nearly sound speed (~ 0.1). Each pair of blue (in-fall) and red (out-moving) 
regions indicates one vortex. Initially there are 5-6 vortices present in the disk since these modes 
have the highest linear growth rates (Paper II). But these vortices have slightly different phase 
speeds going around the disk, which means that vortices will eventually catch up with each other 
and interact. In the end, only one strong vortex left. (A lower resolution version is shown here in 
order to reduce the file size.) 
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Fig. 10. — Angular momentum transport by vortices. Shown are the spatial distribution of a rt f, 
(equation p5|) in the {r, 4>} plane around a vortex from run T5 at time t = 3 (left) and 20 (right) 
orbits. The parameter a r( j, is color-coded. The lower panel shows the azimuthally averaged {a T ,<j>) 
(equation [|l(|) as a function of radius, at the same times as those in the upper panel. Positive 
(a) indicates an outward transport of angular momentum. The transport is peaked at the vortex 
region and remains finite away from the vortex. The strength of the transport still remains high at 
20 orbits. 



Fig. 11. — Similar to Figure 10 except using run Tl at t = 7 (left) and 20 (right) orbits. The 



overall strength of (a) is much smaller than that from T5. 



Fig. 12. — The radial dependence and evolution of (a) (equation jl6)) for runs T1-T8 (top left to 
bottom right). The solid, dotted, dashed and long-dashed curves are at time t = 0, 7, 10, 20 orbits 
for runs T1-T4 and t = 0, 3, 10, 20 orbits for runs T5-T8, respectively. Each row uses the scale 
shown on the left. The angular momentum, on average, is always transported outwards through 
each radius (ring). The transport peaks when the vortices first form but remains steady between 
10 and 20 orbits. 



Fig. 13. — The evolution of the azimuthally averaged physical quantities for run Tl. From top to 
bottom, these are pressure (P), surface density (£), radial mass flux (accretion rate) F%, radial 
angular momentum flux Fj, and angular momentum transport efficiency {a). The solid, dotted 
and dashed curves in each plot are at time t = 0, 10, 20 orbits, respectively. The disk evolution is 
relatively slow, with small changes in averaged quantities (but cf. Figure ^ for changes in azimuthal 
direction) . 



Fig. 14. — Similar to Figure |i~3| but for runs T6 (left) and T8 (right). The solid, dotted, dashed and 
long-dashed curves in each plot are at time t = 0, 3, 10, 20 orbits, respectively. The disk evolution 
is much faster and stronger compared to run Tl. 



Fig. 15. — The growth and saturation of radial velocities for various runs. Runs Tl, Til, T12 and 
T5 have the same initial configurations except for the increasing amplitude A. The saturation level 
increases as the linear growth rate increases, but ceases to do so when the linear growth rate is 
large enough, as shown by runs T5, T6, T7 and T8 (note the difference in their peak values). This 
is because shocks are responsible for the growth saturation in all the high growth rate runs. 
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